NTM analyzed

install packages

#BiocManager::install("clusterProfiler")
#install.packages("enrichplot")
#installed.packages("patchwork")
#BiocManager::install("SingleR")
#BiocManager::install("celldex")
#BiocManager::install("ensembldb")
#BiocManager::install('limma')
#BiocManager::install("DESeq2")
#organism = "org.Mm.eg.db"
#BiocManager::install(organism)
#BiocManager::install("pathview")
#install.packages("Seurat")
#install.packages("remotes")
#remotes::install_github("satijalab/seurat")
#BiocManager::install("rhdf5")
#install.packages("hdf5r")
#BiocManager::install("NeuCA")
#install.packages('devtools')
#install_github("ggjlab/scMCA")
#install.packages("pheatmap")
#BiocManager::install("scater")
#devtools::install_github("sqjin/CellChat")
#BiocManager::install("ComplexHeatmap")
#BiocManager::install("SingleCellExperiment")

Load libraries

In this step, we are loading the libraries which are required for analysis and plotting of the data.

library(Seurat)
library(devtools)
library(ggplot2)
library(patchwork)
library(dplyr)
library(clusterProfiler)
library(enrichplot)
# we use ggplot2 to add x axis labels (ex: ridgeplot)
library(stringr)
library(celldex)
library(ensembldb)
library(SingleR)
library(DESeq2)
library(org.Mm.eg.db)
library(rhdf5)
library(hdf5r)
library(pheatmap)
library(ComplexHeatmap)
library(CellChat)
library(patchwork)
options(stringsAsFactors = FALSE)
library(NMF)
library(ggalluvial)

set the colors of the cell types

my_cols <- c('B cells'='#F8766D','DC'='#39568CFF','Endothelial cells'='#CD9600','Epithelial cells'='#00C19A','Fibroblasts'= "#C77CFF",
    'Macrophages'='#00A9FF','Monocytes'='#0CB702','Stromal cells'='#FF61CC', "ILC" = "grey", "Neutrophils" = "Darkred")

load data for mouse 1 and mouse 2

# Mouse 1
ntm_mouse_1 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/ntm_sample_round_one/outs/"

list.files(ntm_mouse_1)

##  [1] "analysis"                         "cloupe.cloupe"                   
##  [3] "filtered_feature_bc_matrix_ntm_1" "filtered_feature_bc_matrix.h5"   
##  [5] "metrics_summary.csv"              "molecule_info.h5"                
##  [7] "possorted_genome_bam.bam"         "possorted_genome_bam.bam.bai"    
##  [9] "probe_set.csv"                    "raw_feature_bc_matrix"           
## [11] "raw_feature_bc_matrix.h5"         "spatial"                         
## [13] "spatial_enrichment.csv"           "web_summary.html"

ntm_mouse_1_data <- Load10X_Spatial(
  ntm_mouse_1,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "ntm",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)

# Mouse 2
ntm_mouse_2 <- "/Volumes/Seagate_4TB/spatial_ntm_samples/demultiplexed/ntm_sample_round_two/outs/"

list.files(ntm_mouse_2)

##  [1] "analysis"                         "cloupe.cloupe"                   
##  [3] "filtered_feature_bc_matrix_ntm_2" "filtered_feature_bc_matrix.h5"   
##  [5] "metrics_summary.csv"              "molecule_info.h5"                
##  [7] "ntm_only_web_summary.html"        "possorted_genome_bam.bam"        
##  [9] "possorted_genome_bam.bam.bai"     "probe_set.csv"                   
## [11] "raw_feature_bc_matrix"            "raw_feature_bc_matrix.h5"        
## [13] "spatial"                          "spatial_enrichment.csv"

ntm_mouse_2_data <- Load10X_Spatial(
 ntm_mouse_2,
  filename = "filtered_feature_bc_matrix.h5",
  assay = "spatial",
  slice = "ntm",
  filter.matrix = TRUE,
  to.upper = FALSE,
  image = NULL)

data preprocessing

# Mouse 1
ntm_mouse_1_data <- SCTransform(ntm_mouse_1_data, assay = "spatial", verbose = FALSE)

# Mouse 2
ntm_mouse_2_data <- SCTransform(ntm_mouse_2_data, assay = "spatial", verbose = FALSE)

Visualization of lymphoid follicle markers

Mouse 1

p1 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(ntm_mouse_1_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8

Mouse 2

p1 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Cd19", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p2 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Igha", alpha = c(0.1, 1), pt.size.factor = 1.5)
p3 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Bcl6", alpha = c(0.1, 1), pt.size.factor = 1.5)
p4 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Lta", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p5 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Ltb", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p6 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Cd4", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p7 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Mki67", alpha = c(0.1, 1.5), pt.size.factor = 1.5)
p8 <- SpatialFeaturePlot(ntm_mouse_2_data, features = "Cxcr5", alpha = c(0.1, 1.5), pt.size.factor = 1.5)

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8

dimensionality reduction and clustering

# Mouse 1
ntm_mouse_1_analyzed <- RunPCA(ntm_mouse_1_data, assay = "SCT", verbose = FALSE)
ntm_mouse_1_analyzed <- FindNeighbors(ntm_mouse_1_analyzed, reduction = "pca", dims = 1:5)
ntm_mouse_1_analyzed <- FindClusters(ntm_mouse_1_analyzed, verbose = FALSE, resolution = 0.8)
ntm_mouse_1_analyzed <- RunUMAP(ntm_mouse_1_analyzed, reduction = "pca", dims = 1:5)

# Mouse 2
ntm_mouse_2_analyzed <- RunPCA(ntm_mouse_2_data, assay = "SCT", verbose = FALSE)
ntm_mouse_2_analyzed <- FindNeighbors(ntm_mouse_2_analyzed, reduction = "pca", dims = 1:5)
ntm_mouse_2_analyzed <- FindClusters(ntm_mouse_2_analyzed, verbose = FALSE, resolution = 0.8)
ntm_mouse_2_analyzed <- RunUMAP(ntm_mouse_2_analyzed, reduction = "pca", dims = 1:5)

labelling clusters using SingleR pipeline

# First convert the seurat object to single cell experiment object (otherwise the SingleR pipeline will not work)

# Mouse 1
ntm_mouse_1_analyzed.sce <- as.SingleCellExperiment(ntm_mouse_1_analyzed)

# Mouse 2
ntm_mouse_2_analyzed.sce <- as.SingleCellExperiment(ntm_mouse_2_analyzed)

#create reference data

ref.data <- celldex::ImmGenData()

# run singleR pipeline to find the cell types

## Mouse 1
cell_types_mouse_1 <- SingleR(test=ntm_mouse_1_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

## Mouse 2
cell_types_mouse_2 <- SingleR(test=ntm_mouse_2_analyzed.sce, assay.type.test=1, 
    ref=ref.data, labels=ref.data$label.main)

Check if the cell types are labelled with high accuracy

# view cell types

# Mouse 1
table(cell_types_mouse_1$labels)

## 
##           B cells                DC Endothelial cells  Epithelial cells 
##               152                59               868                35 
##       Fibroblasts       Macrophages         Monocytes       Neutrophils 
##               716               476                 7                 1 
##     Stromal cells 
##               601

#checking cell types Mouse 1

plotScoreHeatmap(cell_types_mouse_1)

plotDeltaDistribution(cell_types_mouse_1)

summary(is.na(cell_types_mouse_1$pruned.labels))

##    Mode   FALSE    TRUE 
## logical    2901      14

# Mouse 2
table(cell_types_mouse_2$labels)

## 
##           B cells                DC Endothelial cells  Epithelial cells 
##               124                19               970                28 
##       Fibroblasts       Macrophages         Monocytes       Neutrophils 
##               234               477                 6                 4 
##     Stromal cells 
##               312

#checking cell types Mouse 2

plotScoreHeatmap(cell_types_mouse_2)

plotDeltaDistribution(cell_types_mouse_2)

summary(is.na(cell_types_mouse_2$pruned.labels))

##    Mode   FALSE    TRUE 
## logical    2167       7

merging metadata of cell types to seurat object

# Mouse 1
ntm_mouse_1_analyzed[["ref.data"]] <- cell_types_mouse_1$labels

# Mouse 2
ntm_mouse_2_analyzed[["ref.data"]] <- cell_types_mouse_2$labels

Plot umap and spatial plot to localize each cell type

Mouse 1

p1 <- DimPlot(ntm_mouse_1_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)

p2 <- SpatialDimPlot(ntm_mouse_1_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)

p1 + p2

Mouse 2

p1 <- DimPlot(ntm_mouse_2_analyzed, group.by = c("ref.data"), reduction = "umap",label = F, pt.size = 1.5, label.size = 5, cols = my_cols)

p2 <- SpatialDimPlot(ntm_mouse_2_analyzed, group.by = c("ref.data"), label = FALSE, label.size = 5, cols = my_cols)

p1 + p2

Plot violin plot with top markers associated with B cells and lymphoid follicles

Mouse 1

VlnPlot(ntm_mouse_1_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols =  my_cols)

Mouse 2

VlnPlot(ntm_mouse_2_analyzed, group.by = c("ref.data"), features = c("Cxcr5", "Mki67", "Ltb", "Cd38", "Ccr6", "Tnfrsf13c"), cols =  my_cols)

find markers in each cluster

Mouse 1

B_cell_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_1, n = 10)

##                   p_val avg_log2FC pct.1 pct.2     p_val_adj
## Tnfrsf13c 4.554733e-117   1.255106 0.816 0.142 8.323320e-113
## Ms4a1     3.671873e-111   1.080410 0.724 0.107 6.709980e-107
## Cd79b     2.869009e-105   1.819936 0.914 0.254 5.242828e-101
## Iglc3     9.359354e-105   1.368350 0.868 0.201 1.710328e-100
## Pou2af1   1.805543e-101   1.920921 0.974 0.349  3.299449e-97
## Ptprcap    1.197781e-99   1.093994 0.822 0.170  2.188825e-95
## Cd79a      2.339769e-98   2.275292 0.987 0.425  4.275694e-94
## Ltb        7.343352e-95   1.289891 0.875 0.220  1.341924e-90
## Spib       9.978568e-94   1.299672 0.763 0.157  1.823483e-89
## Cd19       2.272512e-93   1.590420 0.849 0.222  4.152788e-89

write.csv(B_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_mouse_1.csv")

DC_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_1, n = 10)

##                  p_val avg_log2FC pct.1 pct.2    p_val_adj
## H2-Eb1    5.334939e-24  1.0919341 1.000 0.994 9.749067e-20
## C3        4.157713e-23  1.0263028 1.000 0.961 7.597804e-19
## Tmsb4x    1.105303e-21  0.7047462 1.000 1.000 2.019830e-17
## H2-Aa     1.637641e-21  1.0455447 1.000 0.887 2.992626e-17
## Psmb8     4.052254e-21  0.8606870 1.000 0.818 7.405089e-17
## Ccl22     1.634381e-20  0.6154602 0.525 0.123 2.986668e-16
## Inmt      2.315338e-20 -1.8594434 0.458 0.870 4.231048e-16
## Cd52      7.290028e-20  0.9789633 0.983 0.765 1.332180e-15
## H2-Ab1    8.191305e-20  0.8728392 1.000 0.960 1.496879e-15
## Serpina3g 1.069894e-19  1.0868829 1.000 0.829 1.955125e-15

write.csv(DC_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/DC_mouse_1.csv")

Stromal_cells_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_1, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Myh11  6.158379e-36  0.5558271 0.532 0.285 1.125382e-31
## Gsn    3.134160e-33  0.3994196 0.993 0.925 5.727364e-29
## Des    1.691635e-32  0.5752681 0.654 0.432 3.091294e-28
## Ctss   2.689407e-32 -0.6236161 0.952 0.983 4.914621e-28
## Tagln  2.766731e-28  0.4075501 0.461 0.249 5.055924e-24
## Myl9   5.016710e-27  0.4943200 0.496 0.284 9.167536e-23
## Car3   1.236110e-25  1.1447737 0.374 0.189 2.258867e-21
## B2m    3.247562e-25 -0.3333778 0.998 1.000 5.934595e-21
## Tmsb4x 1.472912e-24 -0.3007193 1.000 1.000 2.691599e-20
## Synpo2 1.333461e-23  0.3428623 0.381 0.197 2.436767e-19

write.csv(Stromal_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Stromal_cells_mouse_1.csv")


Endothelial_cells_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cldn5  1.823586e-102  0.6927932 0.992 0.812 3.332422e-98
## Igfbp2  3.821331e-95  0.7524716 0.983 0.778 6.983101e-91
## Epas1   8.978855e-91  0.6125308 0.995 0.899 1.640796e-86
## Ager    1.241455e-86  0.5622126 1.000 0.936 2.268635e-82
## Acer2   8.633030e-86  0.6179743 0.983 0.809 1.577600e-81
## Clic5   2.306495e-85  0.5693432 0.993 0.888 4.214889e-81
## Egfl7   2.202114e-82  0.5948463 0.984 0.834 4.024143e-78
## C3      9.870730e-82 -0.8186092 0.937 0.972 1.803777e-77
## Vegfa   8.421764e-80  0.5778378 0.988 0.865 1.538993e-75
## Ace     3.230574e-78  0.6225789 0.963 0.766 5.903551e-74

write.csv(Endothelial_cells_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Endothelial_cells_mouse_1.csv")

Epithelial_cell_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Gp2     3.057056e-67  1.0931770 0.686 0.044 5.586464e-63
## Tmc5    4.457791e-65  0.7977855 0.600 0.034 8.146167e-61
## Cd177   2.446346e-56  0.8906558 0.629 0.044 4.470452e-52
## B3galt2 1.335568e-51  0.5792686 0.514 0.031 2.440617e-47
## Bpifb1  3.657771e-51  1.9902073 0.800 0.088 6.684210e-47
## Erich3  3.941690e-50  0.9298852 0.800 0.082 7.203045e-46
## Muc5b   1.797909e-47  1.3634905 0.686 0.066 3.285499e-43
## Pklr    4.992440e-47  0.8550121 0.629 0.053 9.123184e-43
## Tff2    1.128400e-43  1.4010201 0.686 0.072 2.062038e-39
## Kcnk12  1.811543e-42  0.7501989 0.657 0.064 3.310413e-38

write.csv(Epithelial_cell_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Epithelial_cell_mouse_1.csv")

Fibroblasts_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_1, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Gsn    2.679271e-39  0.4309985 0.997 0.920 4.896101e-35
## Psap   8.682005e-31 -0.3999956 0.999 1.000 1.586550e-26
## Inmt   1.692121e-30  0.4104036 0.980 0.824 3.092182e-26
## Ccl5   9.335028e-27 -0.6195179 0.538 0.692 1.705883e-22
## Ltbp4  8.209228e-25  0.3411285 0.939 0.819 1.500154e-20
## Ogn    1.929405e-24  0.3335998 0.747 0.543 3.525795e-20
## Ctss   1.099770e-23 -0.5534927 0.973 0.978 2.009719e-19
## Socs1  3.691625e-23 -0.4632449 0.409 0.584 6.746076e-19
## H2-Ab1 2.379039e-22 -0.4181848 0.946 0.965 4.347456e-18
## Fmo2   2.668078e-22  0.3470290 0.957 0.808 4.875645e-18

  write.csv(Fibroblasts_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Fibroblasts_mouse_1.csv")
  
Macrophages_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Macrophages_mouse_1, n = 10)

##                 p_val avg_log2FC pct.1 pct.2     p_val_adj
## Ctss    4.554123e-140  1.2138975 1.000 0.972 8.322204e-136
## Inmt    9.007957e-120 -1.5120459 0.620 0.909 1.646114e-115
## Psap    2.867758e-116  0.8116796 1.000 0.999 5.240540e-112
## B2m     2.043846e-113  0.6627254 1.000 1.000 3.734924e-109
## Gsn     4.765322e-111 -1.0967002 0.836 0.959 8.708150e-107
## Gbp2b   4.006186e-110  1.1110715 0.979 0.754 7.320904e-106
## Tnfaip2 5.875226e-108  1.1302766 0.968 0.771 1.073639e-103
## H2-Eb1  3.943305e-107  0.7819851 1.000 0.993 7.205996e-103
## H2-Ab1  1.595128e-106  0.8231051 1.000 0.953 2.914937e-102
## Cd74    1.150780e-104  0.5097265 1.000 1.000 2.102936e-100

write.csv(Macrophages_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Macrophages_mouse_1.csv")

Monocytes_mouse_1 <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Monocytes", min.pct = 0.25)
  head(Monocytes_mouse_1, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Nlgn3   1.740289e-15  0.3512044 0.286 0.008 3.180203e-11
## Olfr447 6.404532e-15  0.3507122 0.286 0.008 1.170364e-10
## Jph3    4.958409e-13  0.3487454 0.286 0.010 9.060997e-09
## Cma1    6.626389e-12  0.3467812 0.286 0.011 1.210906e-07
## Polr3f  1.413811e-11  0.3462906 0.286 0.011 2.583598e-07
## Zfp846  2.843161e-11  0.3462906 0.286 0.011 5.195592e-07
## Zfp982  3.118175e-11  0.4783202 0.429 0.025 5.698154e-07
## Gm15056 3.191898e-11  0.8252492 0.571 0.047 5.832874e-07
## U2af1l4 8.187911e-11  0.4754203 0.429 0.026 1.496259e-06
## Adgrg3  1.059973e-10  0.3453098 0.286 0.012 1.936995e-06

write.csv(Monocytes_mouse_1, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Monocytes_mouse_1.csv")

Mouse 2

B_cell_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data", ident.1 = "B cells", min.pct = 0.25)
head(B_cell_mouse_2, n = 10)

##                  p_val avg_log2FC pct.1 pct.2    p_val_adj
## Bank1     4.202720e-83  1.3287032 0.726 0.127 5.797232e-79
## Coch      1.488250e-82  0.4608521 0.274 0.008 2.052892e-78
## Tnfrsf13c 1.201551e-81  1.7538640 0.855 0.208 1.657419e-77
## Cxcr5     2.614613e-81  1.0992235 0.653 0.096 3.606597e-77
## Fcrl5     1.683029e-79  0.8158863 0.508 0.051 2.321570e-75
## Ccr6      4.425674e-78  1.5810476 0.790 0.177 6.104775e-74
## Cd79a     1.041844e-77  2.5173171 1.000 0.572 1.437120e-73
## Cd79b     1.195236e-77  2.1139509 0.984 0.383 1.648708e-73
## Cd19      1.761850e-77  1.8569971 0.831 0.207 2.430296e-73
## Spib      2.477439e-76  1.7252838 0.855 0.229 3.417379e-72

write.csv(B_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_mouse_2.csv")

DC_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "DC", min.pct = 0.25)
head(DC_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Xcr1    3.062816e-11  0.3883528 0.368 0.044 4.224848e-07
## Foxd2   8.827878e-11  0.3472199 0.316 0.034 1.217717e-06
## Ccr7    1.753044e-10  0.7015559 0.789 0.206 2.418149e-06
## Ints10  2.099668e-10  0.3006765 0.263 0.025 2.896283e-06
## Ttc39a  3.705338e-10  0.5240484 0.684 0.149 5.111144e-06
## Mettl25 4.552337e-10  0.4161537 0.421 0.063 6.279494e-06
## Cabp4   1.689226e-09  0.3401178 0.316 0.039 2.330119e-05
## Nckap1l 3.071381e-09  0.4610973 0.421 0.070 4.236663e-05
## Adam11  6.130534e-09  0.3530569 0.263 0.030 8.456459e-05
## Cd5     7.717586e-09  0.2941632 0.263 0.030 1.064564e-04

write.csv(DC_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/DC_mouse_2.csv")

Stromal_cells_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Stromal cells", min.pct = 0.25)
head(Stromal_cells_mouse_2, n = 10)

##                 p_val avg_log2FC pct.1 pct.2    p_val_adj
## Myh11    8.950958e-39  0.9616085 0.657 0.331 1.234695e-34
## Tagln    5.600725e-32  0.9005221 0.619 0.328 7.725641e-28
## Acta2    6.307122e-32  1.0824300 0.679 0.385 8.700044e-28
## Ltbp4    6.419828e-30  0.5857921 0.952 0.800 8.855510e-26
## Myl9     1.246700e-27  0.7490902 0.577 0.302 1.719698e-23
## Ptgis    2.324078e-26  0.6232297 0.881 0.669 3.205833e-22
## Eln      8.537455e-25  0.7410893 0.904 0.789 1.177657e-20
## Ppp1r14a 4.089734e-23  0.4813196 0.990 0.826 5.641379e-19
## Gsn      3.191197e-22  0.4398570 0.994 0.899 4.401937e-18
## Ccn2     3.356066e-22  0.6573885 0.853 0.664 4.629358e-18

write.csv(Stromal_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Stromal_cells_mouse_2.csv")


Endothelial_cells_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Endothelial cells", min.pct = 0.25)
head(Endothelial_cells_mouse_2, n = 10)

##                 p_val avg_log2FC pct.1 pct.2     p_val_adj
## Igfbp2  5.232736e-141  1.1102123 0.992 0.727 7.218035e-137
## Cldn5   3.910543e-132  0.9235680 0.996 0.804 5.394203e-128
## Epas1   2.721131e-120  0.8123111 0.999 0.890 3.753528e-116
## Ager    6.818852e-113  0.6825610 1.000 0.950 9.405924e-109
## Pakap.1 1.259013e-109  0.7385255 1.000 0.883 1.736682e-105
## Clic5   1.048858e-105  0.6633032 0.999 0.914 1.446795e-101
## Ace     4.903890e-103  0.7553659 0.982 0.741  6.764425e-99
## C3       1.469702e-98 -0.7505188 0.999 0.997  2.027308e-94
## Tns1     1.850102e-97  0.6977257 0.998 0.880  2.552031e-93
## Egfl7    1.363452e-95  0.6908519 0.995 0.827  1.880745e-91

write.csv(Endothelial_cells_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Endothelial_cells_mouse_2.csv")

Epithelial_cell_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Epithelial cells", min.pct = 0.25)
head(Epithelial_cell_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Erich3  5.496562e-20  0.4546580 0.393 0.039 7.581958e-16
## Cxcl17  2.337558e-17  0.8259062 0.571 0.100 3.224428e-13
## Tex26   2.609954e-17  0.5173289 0.393 0.045 3.600170e-13
## Bpifb1  4.200883e-17  0.5399574 0.286 0.025 5.794698e-13
## Ubxn10  4.527741e-17  0.8378696 0.643 0.127 6.245566e-13
## Lrrc71  9.566381e-17  0.5689175 0.500 0.075 1.319587e-12
## Ankfn1  3.546138e-16  0.3531888 0.321 0.032 4.891543e-12
## Zfp750  7.178780e-16  0.4780801 0.393 0.049 9.902410e-12
## Dynlrb2 1.000800e-15  0.5289973 0.429 0.059 1.380503e-11
## Tmem212 1.184842e-15  0.6631443 0.571 0.103 1.634371e-11

write.csv(Epithelial_cell_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Epithelial_cell_mouse_2.csv")

Fibroblasts_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Fibroblasts", min.pct = 0.25)
head(Fibroblasts_mouse_2, n = 10)

##               p_val avg_log2FC pct.1 pct.2    p_val_adj
## Cst3   3.397128e-24  0.3258673 1.000 1.000 4.685998e-20
## Clec3b 2.544528e-20  0.5677545 0.932 0.746 3.509922e-16
## Gsn    3.803288e-19  0.4445807 0.996 0.903 5.246255e-15
## Psap   3.958205e-17 -0.4074829 1.000 1.000 5.459949e-13
## Ctss   4.568216e-14 -0.4660153 1.000 1.000 6.301398e-10
## Tppp3  5.760751e-14  0.5288847 0.991 0.909 7.946379e-10
## Cybb   1.130674e-13 -0.4894617 0.872 0.926 1.559652e-09
## Lgals3 1.663983e-13 -0.5059102 0.979 0.974 2.295298e-09
## Gbp2b  2.252567e-12 -0.5267407 0.932 0.974 3.107191e-08
## Ctsz   2.496011e-12 -0.4260957 0.927 0.968 3.442998e-08

  write.csv(Fibroblasts_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Fibroblasts_mouse_2.csv")
  
Macrophages_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", min.pct = 0.25)
  head(Macrophages_mouse_2, n = 10)

##                 p_val avg_log2FC pct.1 pct.2     p_val_adj
## Igfbp2  3.409736e-105 -1.5826761 0.627 0.906 4.703390e-101
## Tnfaip2 3.424259e-100  1.1606650 1.000 0.978  4.723423e-96
## Epas1    2.429809e-99 -1.0990461 0.878 0.955  3.351679e-95
## Inmt     3.250042e-99 -1.3557853 0.610 0.906  4.483108e-95
## Gbp2b    3.785450e-99  0.9670132 0.996 0.962  5.221650e-95
## Psap     1.318899e-97  0.6908248 1.000 1.000  1.819289e-93
## Ctss     1.242439e-96  0.8073770 1.000 1.000  1.713820e-92
## Tns1     1.760472e-94 -1.0432070 0.870 0.951  2.428395e-90
## Fmo2     9.261547e-92 -1.1333143 0.618 0.902  1.277538e-87
## Cldn5    8.194109e-91 -1.1740409 0.776 0.922  1.130295e-86

write.csv(Macrophages_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Macrophages_mouse_2.csv")

Monocytes_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Monocytes", min.pct = 0.25)
  head(Monocytes_mouse_2, n = 10)

##                      p_val avg_log2FC pct.1 pct.2    p_val_adj
## Rnf207        7.675060e-54  0.4123782 0.333 0.002 1.058698e-49
## Nup210l       7.992847e-30  0.4090608 0.333 0.004 1.102533e-25
## P2rx2         4.105985e-22  0.4064125 0.333 0.006 5.663795e-18
## Tdrkh         1.039971e-14  0.4011303 0.333 0.010 1.434535e-10
## Ssc5d         3.216632e-12  0.3978387 0.333 0.012 4.437023e-08
## Cdc6          8.956139e-12  0.5442835 0.500 0.028 1.235410e-07
## 2410131K14Rik 4.319930e-10  0.6535795 0.667 0.059 5.958912e-06
## Clstn2        3.541658e-09  0.6435651 0.667 0.065 4.885364e-05
## Crip3         9.169083e-09  0.3899695 0.333 0.018 1.264783e-04
## Isoc2b        5.081043e-08  0.3867032 0.333 0.019 7.008790e-04

write.csv(Monocytes_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Monocytes_mouse_2.csv")

Neutrophils_mouse_2 <- FindMarkers(ntm_mouse_2_analyzed, group.by = "ref.data",  ident.1 = "Neutrophils", min.pct = 0.25)
  head(Neutrophils_mouse_2, n = 10)

##                p_val avg_log2FC pct.1 pct.2    p_val_adj
## Sema5b  4.443551e-25  0.3192712  0.25 0.002 6.129434e-21
## Cyp26a1 4.443551e-25  0.3192712  0.25 0.002 6.129434e-21
## Frmd7   4.443551e-25  0.3192712  0.25 0.002 6.129434e-21
## Zfp469  4.195789e-21  0.3186077  0.25 0.002 5.787671e-17
## Spock1  4.195789e-21  0.3186077  0.25 0.002 5.787671e-17
## Med12   2.364461e-18  0.5809790  0.25 0.003 3.261537e-14
## Cdx2    2.932706e-18  0.3179446  0.25 0.003 4.045375e-14
## Hoxc5   2.932706e-18  0.3179446  0.25 0.003 4.045375e-14
## Cartpt  3.039891e-18  0.3152950  0.25 0.003 4.193226e-14
## Grm5    5.413888e-18  0.5677795  0.50 0.012 7.467918e-14

write.csv(Neutrophils_mouse_2, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/Neutrophils_mouse_2.csv")

B_cell_deseq <- FindMarkers(ntm_mouse_1_analyzed, group.by = "ref.data",  ident.1 = "Macrophages", test.use = "DESeq2")
write.csv(B_cell_deseq, "/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_deseq.csv")

Cell-cell interaction using cell chat

ntm_cellchat <- createCellChat(object = ntm_mouse_2_analyzed, group.by = "ref.data")

## The `data` slot in the default assay is used. The default assay is SCT 
## The `meta.data` slot in the Seurat object is used as cell meta information 
## The cell groups used for CellChat analysis are  B cells DC Endothelial cells Epithelial cells Fibroblasts Macrophages Monocytes Neutrophils Stromal cells

CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)

CellChatDB.use <- CellChatDB
ntm_cellchat@DB <- CellChatDB.use


ntm_cellchat <- subsetData(ntm_cellchat)

## Issue identified!! Please check the official Gene Symbol of the following genes:  
##  H2-BI H2-Ea-ps

ntm_cellchat <- identifyOverExpressedGenes(ntm_cellchat)
ntm_cellchat <- identifyOverExpressedInteractions(ntm_cellchat)


ntm_cellchat <- computeCommunProb(ntm_cellchat)
ntm_cellchat <- filterCommunication(ntm_cellchat, min.cells = 10)

## The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  Monocytes Neutrophils

ntm_cellchat <- computeCommunProbPathway(ntm_cellchat)
ntm_cellchat <- aggregateNet(ntm_cellchat)

groupSize <- as.numeric(table(ntm_cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(ntm_cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(ntm_cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")

# See interaction of each cell type with the another

mat1 <- ntm_cellchat@net$weight
par(mfrow = c(3,4), xpd=TRUE)
for (i in 1:nrow(mat1)) {
  mat2 <- matrix(0, nrow = nrow(mat1), ncol = ncol(mat1), dimnames = dimnames(mat1))
  mat2[i, ] <- mat1[i, ]
  netVisual_circle(mat2, vertex.weight = groupSize, weight.scale = T, edge.weight.max = max(mat1), title.name = rownames(mat1)[i])
}

##Visualize signalling pathways in each cell

# Identify signaling pathways showing significant communications 
ntm_cellchat@netP$pathways

##  [1] "COLLAGEN"    "APP"         "THBS"        "FN1"         "MIF"        
##  [6] "COMPLEMENT"  "GALECTIN"    "MHC-II"      "VEGF"        "LAMININ"    
## [11] "SEMA4"       "ICAM"        "CXCL"        "GAS"         "CCL"        
## [16] "ITGAL-ITGB2" "JAM"         "PECAM1"      "CD22"        "CD45"       
## [21] "PDGF"        "CD52"        "TENASCIN"    "THY1"        "TGFb"       
## [26] "FGF"         "SEMA3"       "CHEMERIN"    "PD-L1"       "LT"         
## [31] "SEMA7"       "SPP1"        "NOTCH"       "GRN"         "SEMA6"      
## [36] "ESAM"        "IL16"        "EPHB"        "PARs"        "CDH"        
## [41] "ncWNT"       "IGF"         "CD39"        "APRIL"       "ICOS"       
## [46] "PDL2"        "NRG"         "XCR"         "IL1"

# Evaluate one of the pathways
pathways.show <- c("MHC-II") 
vertex.receiver = seq(2,5) # a numeric vector. 
?netVisual_aggregate(ntm_cellchat, signaling = pathways.show,  vertex.receiver = vertex.receiver)

# Chord diagram
par(mfrow=c(1,1))
netVisual_aggregate(ntm_cellchat, signaling = pathways.show, layout = "chord")

# Heatmap
par(mfrow=c(1,1))
netVisual_heatmap(ntm_cellchat, signaling = pathways.show, color.heatmap = "Reds")

Compute the contribution of each ligand-receptor pair to the overall signaling pathway

netAnalysis_contribution(ntm_cellchat, signaling = pathways.show)

pairLR.MHC <- extractEnrichedLR(ntm_cellchat, signaling = pathways.show, geneLR.return = FALSE)
LR.show <- pairLR.MHC[1,] # show one ligand-receptor pair

# Hierarchy plot
vertex.receiver = seq(1,4) # a numeric vector
netVisual_individual(ntm_cellchat, signaling = pathways.show,  pairLR.use = LR.show, vertex.receiver = vertex.receiver)

## [[1]]

Visualize cell-cell communication mediated by multiple signaling pathways

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_bubble(ntm_cellchat, sources.use = 1, targets.use = c(2:8), remove.isolate = FALSE)

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# show all the interactions sending from Inflam.FIB
netVisual_chord_gene(ntm_cellchat, sources.use = 1, targets.use = c(2:8), lab.cex = 0.5,legend.pos.y = 30)

Compute the network centrality scores

ntm_cellchat <- netAnalysis_computeCentrality(ntm_cellchat, slot.name = "netP") 

# the slot 'netP' means the inferred intercellular communication network of signaling pathways

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
netAnalysis_signalingRole_network(ntm_cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

ht1 <- netAnalysis_signalingRole_heatmap(ntm_cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(ntm_cellchat, pattern = "incoming")
ht1 + ht2

Identify and visualize communication patterns (outgoing signal)

selectK(ntm_cellchat, pattern = "outgoing")

nPatterns = 3
ntm_cellchat <- identifyCommunicationPatterns(ntm_cellchat, pattern = "outgoing", k = nPatterns)

# river plot
netAnalysis_river(ntm_cellchat, pattern = "outgoing")

# dot plot
netAnalysis_dot(ntm_cellchat, pattern = "outgoing")

Identify and visualize communication patterns (incoming signal)

selectK(ntm_cellchat, pattern = "incoming")

nPatterns = 3
ntm_cellchat <- identifyCommunicationPatterns(ntm_cellchat, pattern = "incoming", k = nPatterns)

# river plot
netAnalysis_river(ntm_cellchat, pattern = "incoming")

# dot plot
netAnalysis_dot(ntm_cellchat, pattern = "incoming")

Gene set enrichment

organism = org.Mm.eg.db

prepare input data for gene set enrichment

# reading in data from deseq2
df = read.csv("/Volumes/Seagate_4TB/spatial_data_analysis/gene_expression/ntm/B_cell_deseq.csv", header=TRUE)

# we want the log2 fold change 
original_gene_list <- df$avg_log2FC

# name the vector
names(original_gene_list) <- df$X

# omit any NA values 
gene_list<-na.omit(original_gene_list)

# sort the list in decreasing order (required for clusterProfiler)
gene_list = sort(gene_list, decreasing = TRUE)

Gene set enrichment

#gene set enrichment function for B cell cluster

gse <- gseGO(geneList=gene_list, 
             ont ="ALL", 
             keyType = "SYMBOL", 
             minGSSize = 3, 
             maxGSSize = 800, 
             pvalueCutoff = 0.05, 
             verbose = TRUE, 
             OrgDb = organism, 
             pAdjustMethod = "none",
             eps = 0.0)

dotplot(gse, showCategory=10, split=".sign") + facet_grid(.~.sign)